R code modified from https://github.com/mjacox/Thermal_Displacement/blob/master/oisst_an.m
First calculate 1982-2011 SST monthly mean. To evaluate SST anomaly:
SST.1 = Monthly SST - SST_30yr_mean
SST.anomaly = SST.1 - Linear trend by monthly series from 1982-202205
Land mask from R/01-2_OISST_land_seaice_mask.Rmd
## The same code from R/01-2_OISST_land_seaice_mask.Rmd(and can check that testing figure)
land_mask <- fread("unzip -p ../data_src/mapfiles/sea_icemask_025d.csv.zip")[, .(lon, lat, x, y, landmask)]
setorder(land_mask, -y, x)
land_mask[,ry:=rowid(x)]
yy <- unique(land_mask[,.(y, ry)]) ## check, NOTE that in OISST.nc file, y is from Northest (89.875) to Southest (-89.875)
landm <- dcast(land_mask, x ~ ry, value.var = "landmask")
xx <- as.numeric(landm[,1]$x)
landm <- as.matrix(landm[,-1])
clim_years = seq(1982,2011) #for climatology
climyrs<- clim_years[length(clim_years)] - clim_years[1] + 1 #30 yrs
monstr <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
Recalculate_30yrs = FALSE
if (Recalculate_30yrs) {
plan(multisession)
options(future.globals.maxSize= 1048576000)
stm <- future_lapply(1:12, function(j) {
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
for (i in clim_years) {
x <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
names(x)[1] <- jstr
x[[1]][which(landm==1)] <- NA_real_
ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
x[[1]][which(ice[[1]]==1)] <- NA_real_
if (i == clim_years[1]) {
stmx <- x
datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
} else {
stmx <- c(stmx, x)
datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
}
}
names(stmx) <- rep(jstr, length(names(stmx)))
stmx <- merge(stmx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time") %>%
aggregate(by=paste0(climyrs, " years"), FUN=mean, na.rm=TRUE)
names(stmx)[1] <- jstr
stmx <- stmx %>% select(jstr) %>% adrop
return (stmx)
})
}
if (Recalculate_30yrs) {
# Note that in an earlier version used "GLDASp4 land mask" but current version use "GLDASp5"
# The differnt version of land mask cause difference in the definition of land areas, so that sea-mask
# and then also cause difference in climatology mean result
save(stm, file="../data_src/stats/sst_1982_2011month_climatology_nc.RData")
} else {
load("../data_src/stats/sst_1982_2011month_climatology_nc.RData")
}
stmlibrary(grid)
plotx <- vector("list", length = 12)
pstrx <- '';
for (j in 1:12) {
plotx[[j]] <- gplotx(stm[[j]], monstr[j], returnx = TRUE, minz = -2.5, maxz = 34.5)
pstrx <- paste0(pstrx, "plotx[[", j, "]],")
}
lay1 <- rbind(c(1,1,2,2,3,3),
c(4,4,5,5,6,6),
c(7,7,8,8,9,9),
c(10,10,11,11,12,12))
evplot <- paste0('grid.arrange(', pstrx, ' layout_matrix=lay1)')gx <- eval(parse(text=evplot))To double check if the plots of climatology are right, use existed data that download from NOAA (NetCDF OISST, 0.5d, monthly mean)
Metadata: https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html
https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc (downloaded 202206 and the data is 1981/12 to 2022/05)
if (!require("ncdf4")) install.packages("ncdf4"); library(ncdf4)
sst_mnfile = "../data_src/stats/sst.mnmean.nc"
if (!file.exists(sst_mnfile)) {
tryCatch({
curl::curl_download("https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc", destfile = sst_mnfile)
}, error = function(e) paste0(e))
}
nx0 <- nc_open(sst_mnfile)
print(nx0) ## sst[lon,lat,time] (Chunking: [360,180,464])
latn1<- ncvar_get(nx0, "lat") # 89.5 - -89.5
lngn1<- ncvar_get(nx0, "lon") # 0 - 359.5
time<- ncvar_get(nx0, "time") # 1981-12.01 - 2020-07-01, we now check every March, August, and December
date<- time %>% as.Date(origin="1800-01-01 00:00:0.0")
maridx <- seq(4, 464, by = 12) #check date[maridx]
augidx <- seq(9, 464, by = 12)
decidx <- seq(13, 464, by = 12)
get_sstmx <- function (midx, minz=-2.5, maxz=34.5) {
dtx <- as.data.table(ncvar_get(nx0, "sst")[,,midx]) %>%
.[,.(sstm=mean(value, na.rm=TRUE)), by=.(V1,V2)] %>%
.[,`:=`(longitude=fifelse(lngn1[V1]>180, lngn1[V1]-360, lngn1[V1]), latitude=latn1[V2])] %>%
.[,.(longitude, latitude, sstm)]
setnames(dtx, 1:2 , c("x", "y"))
gx <- gplotx(dtx, var="sstm", returnx=TRUE, minz=minz, maxz=maxz, no_coord=TRUE, maxlim=180)
gx <- gx +
#gx <- ggplot() +
# geom_tile(data=dt, aes_string(x="longitude", y="latitude", fill="sstm"), alpha=0.8) +
# scale_fill_viridis(limits=c(minsst, maxsst)) +
geom_sf(data = ne_coastline(scale = "large", returnclass = "sf"), color = 'darkgray', size = .3) #+
# coord_sf() +
# xlim(c(-180, 180)) + ylim(c(-90, 90))
return(gx)
}
gx3 <- get_sstmx(maridx)
gx8 <- get_sstmx(augidx)
gx12 <-get_sstmx(decidx)gt3 <- gplotx(as.data.table(stm[[3]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Mar", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)
gt8 <- gplotx(as.data.table(stm[[8]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Aug", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)
gt12 <- gplotx(as.data.table(stm[[12]]) %>% .[,`:=`(x=fifelse(x>180, x-360, x))], "Dec", returnx=TRUE, minz=-2.5, maxz=34.5, maxlim=180)
lay2 <- rbind(c(1,1,2,2),
c(3,3,4,4),
c(5,5,6,6))grid.arrange(gx3, gt3, gx8, gt8, gx12, gt12, layout_matrix=lay2)i <- 2022
j <- 4
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
x <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
names(x)[1] <- "anom"
#dim(x[[1]]) #1440. 720
#dim(landm) #1440, 720
#dim(ice) #1440, 720
x[[1]][which(landm==1)] <- NA_real_
ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
x[[1]][which(ice[[1]]==1)] <- NA_real_
#xt <- matrix()
#length(xt) <-1440 * 720
#dim(xt) <- c(1440, 720)
## if one of x and stm is NA, then the substraction is NA (so numbers of NA values increases)
x[[1]] <- x[[1]] - stm[[j]][[1]]## Just check if plot 2022-04
gx <- gplotx(x, "anom", minz=-3.5, maxz=3.5, return=TRUE) #This anomaly is substrat long-term mean (30yrs, 1982-2011)
y <- read_stars("../data_src/oisst/monthly_anom/202204_anom.nc")
names(y)[1] <- "anom"
gy <- gplotx(y, "anom", minz=-3.5, maxz=3.5, return=TRUE)
layt <- rbind(c(1,1,2,2))
grid.arrange(gx, gy, layout_matrix=layt)yrng <- seq(1982,2022)
endt <- "2022-05-22"
trackdate <- seq.Date(as.IDate(endt)-6, as.IDate(endt), by="day")
curryr <- year(as.IDate(endt))
currmo <- month(as.IDate(endt))
NEED_RETEST_PRED <- FALSE
if (NEED_RETEST_PRED) {
j=4 #just testing, arbitrary select a month
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
for (i in yrng) {
if (i==curryr & j>(currmo-1)) break
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
if (i==yrng[1]) {
stdrx <- z
datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
} else {
stdrx <- c(stdrx, z)
datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
}
}
names(stdrx) <- rep(jstr, length(names(stdrx)))
xt <- merge(stdrx)
st_set_dimensions(xt, 3, values = as.POSIXct(datex), names = "time")
seqx <- rbind(c(`377`= as.numeric(xt[[1]][377,41,])), ## arbitrary select a point, check which(!is.na(xt[[1]][,41,]))
c(`1112`= as.numeric(xt[[1]][1112,41,])),
c(`test`= sort(rnorm(41)) + rnorm(41)))
colnames(seqx) <- gsub("-", "_", datex)
rownames(seqx) = c("0377_41", "1222_41", "test")
yhat <- stats::predict(lm(seqx[3,] ~ seq_along(datex)))
seqx[1, 28:30] <- NA_real_ #insert NA
seqx[2, 32:34] <- NA_real_ #insert NA
predy <- function(x) { return (stats::predict(lm(x ~ seq_along(datex)))) }
haty <- apply(seqx, MARGIN=1, FUN=function(y){ #predy use predict will auto exclude NA, so that out length is shorter
if (all(is.na(y))) return(y)
yh <- as.numeric(predy(y))
y[!is.na(y)] <- yh
return(c(as.numeric(y)))
}) %>% as.data.frame()
all.equal(haty$test, as.numeric(yhat)) #TRUE
#test using package pracma detrend function
if (!require("pracma")) install.packages("pracma"); library(pracma)
ydet <- as.numeric(detrend(seqx[3,], 'linear'))
all.equal(ydet, as.numeric(seqx[3,]) - as.numeric(yhat)) #TRUE
#test using package pracma detrend function
if (!require("pracma")) install.packages("pracma"); library(pracma)
plot_pred <- function(x, y, yhat=NULL) {
if (is.null(yhat)) {
yhat <- stats::predict(lm(y ~ x))
} else {
yh <- stats::predict(lm(y ~ x))
if (!isTRUE(all.equal(as.numeric(yh), yhat))) {
print(paste0("Not equal yhat with lm: ", paste(yh, collapse=",")))
}
}
plot(x=x, y=y)
abline(lm(y ~ x))
points(x=x, y=yhat, col="red", pch=4)
points(x=x, y=c(y-yhat), col="green", pch=13)
}
plot_pred(seq_along(datex), seqx[3,], haty$`test`)
}Recalculate_anom <- FALSE
if (Recalculate_anom) {
for (i in yrng) {
mmx <- fifelse(i==curryr, currmo-1L, 12L)
for (j in 1:mmx) {
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
filet <- paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc")
if (!file.exists(filet)) {
#z<- read_stars(paste0("../data_src/oisst/monthly_anom/", i, monj, "_anom.nc")) #Old vers use anom, not sst, so no subtract 30yr mean
z <- read_stars(paste0("../data_src/oisst/monthly_sst/", i, monj, "_sst.nc"))
names(z)[1] <- "anom"
z[[1]][which(landm==1)] <- NA_real_
ice <- read_stars(paste0("../data_src/oisst/monthly_icemask/", i, monj, "_icemask.nc"))
z[[1]][which(ice[[1]]==1)] <- NA_real_
z[[1]] <- z[[1]] - stm[[j]][[1]] #Old-version without this substraction while using _anom directly
write_stars(z, filet)
}
}
}
}## Just check if plot 2022-04
if (Recalculate_anom) {
gx2 <- gplotx(z, "anom", minz=-3.5, maxz=3.5, return=TRUE) #This anomaly is substrat long-term mean (30yrs, 1982-2011)
#y <- read_stars("../data_src/oisst/monthly_anom/202204_anom.nc")
#names(y)[1] <- "anom"
#gy <- gplotx(y, "anom", minz=-3.5, maxz=3.5, return=TRUE)
layt <- rbind(c(1,1,2,2))
grid.arrange(gx2, gy, layout_matrix=layt)
}## It's hard/slow to detrend each x,y along almost 30yr (30x12) year trend. Must read 1440x720x30x12 at once...
Recalculate_anom <- FALSE
yrs <- yrng[length(yrng)] - yrng[1] + 1
if (Recalculate_anom) {
for (j in 1:12) {
monj <- fifelse(j<10, paste0("0",j), paste0(j))
jstr <- monstr[j]
for (i in yrng) {
if (i==curryr & j>(currmo-1)) break
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", i, monj, "_anom.nc"))
if (i==yrng[1]) {
stdrx <- z
datex<- as.Date(paste0(i, monj, "01"), format="%Y%m%d")
} else {
stdrx <- c(stdrx, z)
datex<- c(datex, as.Date(paste0(i, monj, "01"), format="%Y%m%d"))
}
}
names(stdrx) <- rep(jstr, length(names(stdrx)))
print(paste0("Now in Year-month: ", i," - ", monj, " and have length stdrx: ", length(datex)))
predy <- function(x) { return (stats::predict(lm(x ~ seq_along(datex)))) }
trd <- merge(stdrx) %>% st_set_dimensions(3, values = as.POSIXct(datex), names = "time") %>%
aggregate(by=paste0(yrs, " years"), FUN = function(y) {
if (all(is.na(y))) return(list(as.numeric(y)))
y[!is.na(y)] <- as.numeric(predy(y)) #Seems if return value has equal length, the result will be simplified by aggregate if using as.matrix
return(list(as.numeric(y))) #Old-version bug: not predy(y), should be y. See previous testing chunk
})
tk <- array(unlist(trd[[1]]), dim = c(length(datex),1440,720)) #bug fix: NOT each predict have the same length result (if input has NA)
print(paste0("Predict ok with dim(tk): ", paste0(dim(tk), collapse=",")))
for (k in seq_along(yrng)) {
if (yrng[k]==curryr & j>(currmo-1)) break
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", yrng[k], monj, "_anom.nc"))
#notna1 <- which(!is.na(z[1][[1]]))
#notna2 <- which(!is.na(tk[k,,]))
#notna <- intersect(notna1, notna2)
#z[[1]][notna] <- z[[1]][notna]- tk[1,,][notna]
z[[1]] <- z[[1]]- tk[1,,] #Any of z[[1]] or tk[1,,] NA will cause NA here, it makes sense.
names(z)[1] <- "anom"
filet <- paste0("../data_src/oisst/monthly_anom_detrend/", yrng[k], monj, "_anom.nc")
write_stars(z, filet)
print(paste0("Now write Year-month: ", yrng[k], " - ", monj, " ok"))
}
}
} yk <- read_stars("../data_src/oisst/monthly_anom/202006_anom.nc")
#zk <- read_stars(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "06", "_anom.nc"))
z <- read_stars(paste0("../data_src/oisst/monthly_anom_icemask/", 2020, "06", "_anom.nc"))
names(yk)[1] <- "anom"
#names(zk)[1] <- "anom"
names(z)[1] <- "anom"
print(range(na.omit(as.data.table(yk)$anom))) #-5.911333 10.598333
#print(range(na.omit(as.data.table(zk)$anom))) #-5.39709 16.39870 #old -8.227043 11.727673
print(range(na.omit(as.data.table(z)$anom))) #-7.724047 10.483911 #old -8.361936 6.181935
gyk<- gplotx(yk,"anom", minz = -6, maxz = 11.0, returnx=TRUE)
gzk<- plot_anom(paste0("../data_src/oisst/monthly_anom_detrend/", 2020, "06", "_anom.nc"), #zk,
"anom", #minz = -8, maxz = 11.0,
returnx=TRUE)
gz <- gplotx(z, "anom", minz = -6, maxz = 16.5, returnx=TRUE) layt <- rbind(c(1,1),c(2,2),c(3,3))
grid.arrange(gyk, gzk, gz, layout_matrix=layt)#gyk
#gzk
#gzif (!require("layer")) {
if (!require("remotes")) install.packages("remotes"); library(remotes)
remotes::install_github("marcosci/layer")
library(layer)
}
Reload_sst_nc <- FALSE
if (Reload_sst_nc) {
sst1 <- read_stars(paste0("../data_src/oisst/monthly_sst/202002_sst.nc"))
names(sst1)[1] <- "sst"
sst202002 <- as.data.table(sst1) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
#.[,.(longitude, latitude, heatwave)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
sst2 <- read_stars(paste0("../data_src/oisst/monthly_sst/201002_sst.nc"))
names(sst2)[1] <- "sst"
sst201002 <- as.data.table(sst2) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
sst3 <- read_stars(paste0("../data_src/oisst/monthly_sst/200002_sst.nc"))
names(sst3)[1] <- "sst"
sst200002 <- as.data.table(sst3) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
sst4 <- read_stars(paste0("../data_src/oisst/monthly_sst/199002_sst.nc"))
names(sst4)[1] <- "sst"
sst199002 <- as.data.table(sst4) %>%
.[,`:=`(longitude=fifelse(x>180, x-360, x), latitude=y)] %>%
.[x>=100 & x<=180 & latitude>=-45 & latitude <= 35,] %>% #x<=300
.[,`:=`(x=NULL, y=NULL)] %>%
st_as_sf(coords=c("longitude", "latitude"))
save(sst199002, sst200002, sst201002, sst202002, file="../data_src/stats/sst_1990_2020_stack.RData")
} else {
load("../data_src/stats/sst_1990_2020_stack.RData")
}
tl0 <- layer::tilt_map(sst199002)
tl1 <- layer::tilt_map(sst200002, y_shift=75)
tl2 <- layer::tilt_map(sst201002, y_shift=150)
tl3 <- layer::tilt_map(sst202002, y_shift=225)
map_list <- list(tl0, tl1, tl2, tl3)#palettem: https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html
layer::plot_tiltedmaps(map_list, palette = "turbo")